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Abstract 

A  probabilistic  analysis  of  the  fatigue  life  of  specimens  subject  to  fretting  fatigue  was 
developed.  A  mechanics  based  fretting  life  analysis  was  applied  that  employed  the  local  stress 
gradient  at  the  edge  of  contact.  The  random  variables  in  the  analysis  included  the  initial  crack 
size,  coefficient  of  friction,  crack  growth  rate  law,  and  the  contact  pad  profile.  The  variation  in 
pad  profiles  was  determined  through  measurement  of  seventy-seven  machined  pads. 
Distributions  for  the  other  random  variables  were  obtained  using  previously  generated  test  data. 
A  probabilistic  fatigue  analysis  was  applied  using  Monte  Carlo  sampling  to  determine  the 
statistics  (mean  and  standard  deviation)  of  the  fatigue  life  prediction  and  was  compared  to 
fretting  fatigue  test  data.  Several  qualitative  and  quantitative  sensitivity  methods  were  applied  to 
the  results  including  the  calculation  of  the  probabilistic  sensitivities  (partial  derivatives  of  the 
fatigue  life  statistics  with  respect  to  the  input  probability  density  function  parameters)  via  linear 
regression  and  finite  difference. 


1  Materials  Research  Engineer,  2230  Tenth  St,  Ste  1,  Wright-Patterson  AFB,  OH  45433-7817,  USA,  Phone:  937- 
255-5438,  FAX:  937-656-4840,  E-mail:  patrick.golden@wpafb.af.mil 


1 


Keywords 

Probabilistic  Analysis,  Sensitivities,  Fretting  Fatigue,  Crack  Propagation 

1.  Introduction 

Fretting  is  a  problem  in  many  aerospace  applications  including  the  blade  to  disk 
attachment  in  turbine  engines.  Two  fretting  modes  often  contribute  to  damage  in  fretting:  gross 
slip  when  the  two  surfaces  slide  resulting  in  wear,  and  partial  slip  when  the  two  surfaces  are 
nominally  stuck  together  except  for  a  small  slip  zone  at  the  edge  of  contact.  The  surface  damage 
and  wear  caused  by  fretting  is  a  costly  maintenance  burden  and  when  combined  with  the  very 
high  local  contact  stresses  due  to  fretting,  it  can  result  in  disk  or  blade  cracking  and  the  potential 
for  catastrophic  failure.  Understanding  of  these  damage  mechanisms  associated  with  fretting  and 
accurate  life  prediction  tools  is  critical  to  the  sustainment  of  aerospace  components. 

Many  approaches  have  been  employed  to  aid  in  the  life  prediction  of  components  with 
fretting.  One  approach  is  the  use  of  material  stress  or  strain-life  curves  with  a  “knockdown 
factor”  determined  through  testing  of  specimens  with  fretting  behavior  similar  to  a  component, 
or  with  actual  or  simulated  components.  This  has  the  disadvantage,  however,  of  being  applicable 
to  only  a  specific  problem.  Other  approaches  include  development  of  critical  lifing  parameters 
that  account  for  combinations  of  slip,  applied  stress,  and/or  strain  [1,2].  These  models  are  more 
general,  but  are  typically  calibrated  to  specimen  data  which  may  not  fully  account  for  the 
different  geometry  and  stress  gradients  present  in  a  component.  Others  have  applied  an  approach 
that  uses  the  local  contact  stress  calculation  combined  with  traditional  multiaxial  lifing  tools, 
such  as  the  Smith-Watson-Topper  or  equivalent  stress  parameters  [3,4],  and  often  combined  with 
fracture  mechanics  fatigue  crack  growth  [5-7].  This  approach  is  well  suited  to  fretting  fatigue 
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problems  that  are  dominated  by  partial-slip  fretting  since  the  surface  damage  and  the  high  crack 
growth  driving  stresses  are  both  localized  at  the  edge  of  contact. 

A  probabilistic  analysis  is  often  applied  to  problems  to  help  determine  the  effect  of 
variability  of  the  model  input  parameters  on  the  model  outputs.  Prior  work  on  statistical  or 
probabilistic  analysis  in  fretting  includes  modeling  of  the  variability  in  the  contact  surface  profile 
by  Kumari  and  Farris  [8].  Here,  the  measured  profiles  of  the  indentors  in  fretting  fatigue  tests 
were  measured  and  statistically  described  and  carried  through  the  life  prediction  models  to 
estimate  the  expected  variability  in  stress  and  life.  Other  probabilistic  fretting  analyses  include 
work  on  fretting  fatigue  of  riveted  lap  joints  [9]  and  on  rolling  contact  [10].  In  any  probabilistic 
analysis  it  is  crucial  to  understand  which  are  the  important  random  variables  driving  the  results 
of  the  analysis.  Importance  is  a  combination  of  those  that  to  which  the  results  are  the  most 
sensitive  and  that  are  alterable  with  the  least  cost. 

The  objective  of  the  current  work  was  to  develop  and  demonstrate  a  probabilistic  fretting 
fatigue  lifing  approach  for  a  dovetail  fretting  experimental  configuration  that  included  a  broad 
range  of  input  random  variables,  and  to  apply  efficient  sensitivity  methods  to  determine  the 
relative  importance  of  the  input  variables.  A  previously  demonstrated  stress  fracture  mechanics 
based  model  was  chosen  for  life  prediction  from  small  fretting  crack  to  failure  was  adapted  to 
this  effort.  Input  probability  density  functions  (PDF’s)  were  developed  using  available 
laboratory  data.  Measures  were  applied  to  the  analysis  to  determine  qualitatively  and 
quantitatively  the  random  variable  parameters  that  resulted  in  the  highest  results  sensitivity. 

2.  Deterministic  Analysis 

A  series  of  fretting  experiments  was  previously  conducted  [11]  to  improve  understanding 
of  fretting  behavior  in  Ti-6A1-4V  and  to  assess  life  prediction  models.  The  geometry  of  the 


3 


fretting  samples  was  a  dovetail  shaped  specimen  that  was  designed  to  represent  the  attachment 
between  a  turbine  engine  blade  and  disk.  The  tests  were  conducted  at  room  temperature,  which 
is  consistent  with  the  operating  conditions  of  a  fan  disk.  The  contact  interface  was  bare  Ti-6A1- 
4V  on  bare  Ti-6A1-4V,  but  several  coatings  and  residual  stress  surface  treatments  applied  to  the 
contact  interface  were  also  tested.  The  analysis  in  the  current  work  was  limited  to  the  bare  Ti- 
6A1-4V  tests.  A  schematic  of  the  fretting  fatigue  test  rig  modeled  in  this  analysis  is  shown  in 
Figure  1.  The  experimental  setup  consists  of  the  dovetail  specimen  (one-half  of  the  specimen  is 
shown  in  the  schematic),  two  fretting  pads,  and  a  steel  fixture.  Here,  Ej  and  V/  are  the  Young’s 
modulus  and  Poissons  Ratio  of  the  specimen,  respectively,  and  E2  and  v2  aret  eh  Young’s 
Modulus  and  Poisson’s  ratio  of  the  pad,  respectively.  Both  the  dovetail  specimen  and  pads  were 
machined  from  Ti-6A1-4V  with  a  thickness  of  7.62  mm,  Young’s  modulus  £=116  GPa,  and 
Poisson’s  ratio  v=  0.31.  The  fretting  pads  are  held  in  the  steel  fixture  at  a  45°  flank  angle,  and 
the  dovetail  specimen  is  pulled  with  a  cyclic  load,  F,  into  the  fretting  pads.  Both  the  dovetail 
specimen  and  pads  were  machined  from  Ti-6A1-4V  with  a  thickness  of  7.62  mm,  modulus  E  = 
116  GPa,  and  Poisson’s  ratio  v=  0.31.  The  nominal  fretting  pad  geometry  was  a  3  mm  flat  with 
3  mm  blending  radii.  The  normal,  P ,  and  shear,  Q ,  contact  forces  were  measured  indirectly  from 
strain  gage  instrumentation  of  the  experiments.  More  details  of  the  experimental  setup  and 
instrumentation  can  be  found  in  Golden  and  Nicholas  [11].  10  tests  were  run  at  values  of  Fmax 
ranging  from  18  kN  to  30  kN  at  a  load  ratio  R  =  0.1. 
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Experiment  Analysis 


Figure  1:  Schematic  of  the  dovetail  experiment  and  the  analysis. 


A  fracture  mechanics  based  fretting  life  prediction  analysis  for  dovetail  specimens  was 
developed  and  demonstrated  using  these  and  other  experimental  test  data  and  is  described  in 
Golden  and  Calcaterra  [12].  The  objective  of  this  prior  work  was  to  demonstrate  lifing  methods 
that  could  be  applied  to  more  complex  3-D  structures  found  in  turbine  engines,  however,  the 
current  problem  can  be  simplified  to  a  2-D  problem.  The  analysis  was  broken  into  two  parts:  a 
finite  element  method  (FEM)  analysis  and  a  2-D  numerical  contact  stress  analysis.  A  finite 
element  analysis  is  needed  in  this  problem  to  determine  two  sets  of  quantities.  The  first  quantity, 
described  here,  relates  to  the  behavior  of  the  contact  force  history.  The  second,  which  can  be 
obtained  from  the  same  FEM  model,  is  the  bulk  stress  as  shown  in  Figure  1  and  is  described 
below.  The  dovetail  geometry  was  modeled  using  a  nonlinear  contact  FEM  model.  This 
analysis  yielded  the  contact  force  history  for  P  and  Q  and  an  example  is  plotted  in  Figure  2.  This 
example  shows  a  single  loading  and  unloading  starting  from  zero  load.  As  the  load  was  applied, 
the  contact  was  initially  in  sliding  (gross  slip)  and  the  contact  forces  followed  the  dashed  line 
defined  by  the  equation 
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Q  =  liP 


(1) 


where  //  is  the  coefficient  of  friction.  Upon  load  reversal,  the  contact  load  path  changes 
slope  to  the  partial  slip  slope,  a,  which  was  positive  in  this  case.  In  real  components  the  slope, 
a ,  could  range  from  a  negative  to  a  positive  value  including  an  infinite  or  vertical  slope  and  is 
dependent  on  the  component  geometry  and  compliance.  The  contact  forces  will  follow  the  slope 
a  for  both  increasing  and  decreasing  loads  according  to  the  equation 

A  Q  =  aAP  (2) 

until  the  dashed  lines  defining  friction  are  reached  according  to  Equation  1  for  increasing 
load  or  the  equation 

Q  =  —pP  (3) 

for  decreasing  load.  Therefore,  the  two  key  inputs  needed  for  prediction  of  the  contact 
force  history  for  a  given  applied  load  history  were  ju  and  a.  The  coefficient  of  friction,  //,  is  a 
property  that  can  only  be  measured,  however,  a  was  predicted  from  the  FEM  analysis  and 
confirmed  by  experimental  measurement.  Once  these  2  parameters  were  known,  the  contact 
forces  could  be  determined  for  a  given  applied  force  without  additional  FEM  analysis  through  a 
procedure  developed  by  Gean  and  Farris  [13].  The  procedure  is  outlined  as  follows  and  has  also 
been  used  on  a  dovetail  application  in  [14].  The  remote  applied  force,  F,  and  contact  forces  P 
and  Q  must  satisfy  the  equation 

F  =  P  sin  9  +  Q  cos  6  (4) 

where  0  is  the  dovetail  flank  angle.  Additionally,  the  only  valid  values  of  P  and  Q  were 
bounded  by  Equations  1  and  X.  If  the  applied  force  F  is  increasing,  the  Q  versus  i  path  will 
always  move  up  with  either  a  slope  //  or  a:  n  if  the  current  location  is  already  on  the  line  defined 
by  Equation  1,  and  a  if  the  current  location  is  below  that  line.  Likewise,  if  F  is  decreasing,  the  Q 
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versus  P  path  will  always  move  to  the  left  with  slope  -//if  the  current  location  is  already  on  the 
line  defined  by  Equation  3,  or  down  with  a  slope  a  if  the  current  location  is  above  that  line. 


P  (N/mm) 

Figure  2:  A  typical  plot  of  Q  versus  P  showing  the  partial  slip  slope  a  and  the  boundaries 

defined  by  // 

Once  the  contact  force  history  was  known,  the  contact  stresses  were  calculated  using  a  2- 
D  numerical  contact  stress  analysis.  Figure  1  shows  a  schematic  of  the  equivalent  contact 
geometry,  defined  by  the  gap  function,  h(x),  used  in  this  analysis.  In  the  experiments  described 
here,  h(x )  was  simply  the  profile  of  the  fretting  pads  that  are  pressed  into  contact  with  the  flat 
dovetail  specimens.  Rather  than  use  the  prescribed  profile  in  the  analysis,  the  as-machined 
profiles  were  measured  using  a  contacting  profilometer.  The  analytical  tool,  CAPRI  (Contact 
Analysis  for  Profiles  of  Random  Indenters)  described  in  McVeigh  et  al.  [15],  was  developed  to 
solve  the  singular  integral  equation  that  defines  this  contact  problem  using  any  reasonable 
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function,  h(x),  as  measured  from  the  profilometry.  The  output  of  CAPRI  is  the  pressure,  p(x), 
and  shear,  q(x),  tractions  due  to  the  applied  contact  forces.  An  example  of  these  tractions  is 
shown  in  Figure  3.  Next,  CAPRI  was  used  to  determine  the  subsurface  stress  field,  Oij(x,y), 
under  the  contact  due  to  p(x)  and  q(x).  Note  the  peaks  located  near  the  edges  of  contact  in  Figure 
3.  These  peaks  in  the  pressure,  shear,  and  the  underlying  subsurface  stresses  were  driven  by  the 
geometry  near  the  edges  of  contact  and  typically  result  in  crack  initiation  and  growth  from  the 
edge  of  contact. 


Figure  3:  Pressure  and  shear  tractions  from  CAPRI  at  maximum  applied  load 

The  last  step  in  the  analysis  was  the  fracture  mechanics  based  fatigue  crack  growth 
prediction.  Here,  the  total  stresses  along  the  expected  crack  growth  path  were  used  to  calculate 
the  stress  intensity  factor  range,  A K.  To  get  the  total  stress  solution  both  the  contact  stresses 
from  CAPRI  and  the  bulk  stress  as  shown  in  Figure  1  were  superposed.  The  bulk  stress  was  the 
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stress  distribution  induced  in  the  component  due  to  the  remotely  applied  loading  and  the 
geometry  of  the  component.  A  procedure  to  obtain  the  bulk  stress  distribution  for  a  dovetail 
geometry  was  described  by  Golden  and  Calcaterra  [12]  which  necessitated  both  the  FEM  model 
and  CAPRI  results.  Once  the  full  subsurface  stress  distribution  was  known,  the  mode  I  stress 
distribution  could  be  extracted  along  the  expected  crack  path  at  the  edge  of  contact.  Since  this 
stress  distribution  was  nonlinear,  weight  function  stress  intensity  factor  solutions  were  applied 
[16].  The  crack  propagation  life  was  then  integrated  using  a  crack  growth  model  of  the  form 

—  =  /( AK)  (2) 

dN  ’ 

where  da/dN  is  the  crack  growth  rate  and  the  crack  growth  model/is  discussed  later. 
Integration  of  life  was  performed  using  the  Euler  method  with  small  step  size  starting  from  a 
small  initial  fretting  crack  typically  25  mm  in  depth,  until  fracture.  Crack  growth  properties 
from  previous  testing  on  the  same  Ti-6A1-4V  material  used  in  the  US  Air  Force  High  Cycle 
Fatigue  program  [17]  were  applied. 

3.  Statistical  Analysis 

The  random  variables  for  the  probabilistic  fretting  analysis  performed  in  this  study  were 
chosen  from  the  key  input  variables  of  the  deterministic  analysis  described  above.  These  key 
input  variables  included  the  initial  crack  size,  the  coefficient  of  friction,  the  slope  a  that  defines 
the  contact  forces  during  partial  slip,  the  crack  growth  law  parameters,  and  the  parameters 
defining  the  shape  of  the  pad  profile.  PDF’s  were  developed  for  each  of  these  random  variables 
from  measurements  and  test  data  available  both  from  the  dovetail  fretting  experiments  and  from 
other  fatigue  and  fatigue  crack  growth  experiments.  The  details  of  the  development  of  these 
PDF’s  are  described  below.  A  summary  of  the  resulting  PDF’s  are  shown  in  Table  1. 
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Initial  crack  size 


The  deterministic  fretting  fatigue  life  prediction  analysis  described  above  was  a  fracture 
mechanics  based  crack  propagation  calculation.  To  perform  this  analysis  an  initial  crack  size 
was  required.  Previous  work  [16]  has  shown  that  the  initial  fretting  cracks  tend  to  be  shallow, 
low  aspect  ratio  (ale)  surface  cracks  that  often  appear  to  be  several  micro-cracks  that  have 
coalesced.  Here,  a  is  the  surface  crack  depth  and  c  is  the  half  surface  length.  Investigations  of 
naturally  initiated  cracks  in  smooth  bar  fatigue  specimens  in  Ti-6A1-4V  with  the  same 
microstructure  [18]  have  revealed  that  the  cracks  appear  to  initiate  almost  exclusively  at  primary 
alpha  grains  at  the  surface  of  the  specimen.  These  naturally  initiated  cracks  in  the  primary  alpha 
grains  form  readily  identifiable  facets  on  the  specimen  fracture  surface.  Measurements  of  these 
initial  cracks  were  made  in  a  prior  study  [18]  on  fatigue  variability  of  Ti-6A1-4V,  where  repeated 
tests  were  conducted  at  a  few  constant  amplitude  stress  levels.  The  resulting  distribution  of 
primary  alpha  grain  facets  was  used  here  to  represent  the  expected  variability  in  initial  fretting 
crack  depth.  It  is  therefore  assumed  in  this  analysis  that  the  naturally  initiated  fretting  cracks 
also  form  at  the  surface  primary  alpha  grains,  albeit  at  a  different  local  stress  state  and  surface 
condition  as  in  the  prior  study.  Since  it  has  been  shown  that  typically  multiple  surface  cracks 
form  and  coalesce  in  fretting,  a  constant,  low  aspect  ratio  of  ale  =  0.2  was  used.  The  mean  and 
standard  deviation  of  the  measured  naturally  initiated  crack  sizes  is  listed  in  Table  1.  It  was 
determined  that  a  lognormal  distribution  was  the  best  fit  to  the  data. 

Coefficient  of  friction  and  the  partied  slip  slope 

As  described  above  and  depicted  in  Figure  2,  the  coefficient  of  friction  and  the  partial  slip 
slope  of  Q  versus  P  were  measured  during  each  test.  The  accuracy  of  these  measurements, 
however,  was  limited  to  the  accuracy  of  the  strain  gage  instrumentation  of  the  tests  and  the  finite 
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element  analysis  that  was  used  to  convert  fixture  strain  to  applied  load.  Additionally,  variability 
in  friction  was  expected  from  test  to  test.  Both  the  measurement  uncertainty  and  inherent 
variability  in  the  tests  needed  to  be  captured  in  the  probabilistic  model.  To  achieve  this,  data  was 
collected  from  twenty-three  fretting  tests  all  of  which  were  conducted  under  similar  loading 
conditions.  All  tests  were  bare  Ti-6A1-4V  on  Ti-6A1-4V  contact  (no  coatings)  and  all  were 
loaded  at  R  =  0.1.  The  measured  values  of  friction  coefficient  and  partial  slip  slope,  a,  were 
found  to  be  correlated.  The  uncertainty  in  the  expected  values  of  friction  and  partial  slip  slope 
along  with  their  correlation  was  modeled  using  a  multivariate  normal  distribution.  In  order  to 
avoid  the  possibility  of  an  infinite  slope,  a,  the  value  1/a  is  used  in  the  analysis  procedure  and 
tabulated  in  Table  1. 

Fatigue  crack  growth  parameters 

The  crack  growth  rate  curves  used  for  crack  propagation  predictions  were  based  on 
fatigue  crack  growth  tests  previously  conducted  [17].  Four  tests  conducted  at  R  =  0.1  and  -1 
from  two  different  labs  were  fit  to  a  bilinear  crack  growth  model,  where,  da/dN  is  the  crack 
growth  rate  and  A K  is  the  stress  intensity  factor  range.  These  data  are  plotted  in  Figure  4.  The 
Walker  model  [19]  was  used  to  collapse  the  data  from  the  different  R  values  using  an  equation  of 
the  form 

da/dN  =  Ct  [AK(1  -  R)(™-D]ni  A  K  <b 
da/dN  —  C2[AK(1  —  R)(m-1)]n2  AK  >  b 

where  b  is  the  intersection  point  for  regions  1  and  2,  defined  as 

^  _  l°§io  Cj  ~~  log  ip  C2  ^ 

n2  -  nx 

and  the  values  of  the  Walker  exponent  m  were  determined  in  the  previous  program  [9]. 
Different  values  were  used  for  R  <  0,  m  =  0.275,  and  R  >  0,  m  =  0.72.  The  statistics  of  the  curve 
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da/dN  (m/cycle) 


fit  were  modeled  using  a  multivariate  normal.  The  coefficients  C,  and  the  exponents  n,  were 


found  to  be  correlated  for  each  segment  of  the  bilinear  crack  growth  model.  The  mean,  standard 
deviation,  and  correlation  coefficients  are  summarized  in  Table  1.  The  values  of  the  model 
parameters  resulted  units  of  MPaVm  for  A K,  and  m/cycle  for  da/dN.  The  da/dN  versus  A K 
curves  for  one-hundred  randomly  generated  realizations  of  C,  and  /?,  were  also  plotted  in  Figure 
4. 


AK  ff  (MPa-m1/2) 

Figure  XX:  Plot  of  da/dN  versus  A K  data  and  100  realizations  of  the  bilinear  model. 
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Pad  profile 

The  nominal  geometry  for  the  contact  pad  is  the  thick  black  line  shown  in  Figure  5.  The 
prescribed  profile  is  the  geometry  that  was  requested  on  the  part  drawings  provided  to  the 
machine  shop.  The  central  segment  of  the  prescribed  profile  is  flat  with  a  length  of  3.00  mm. 

The  flat  section  is  flanked  by  blending  radii  of  3.0  mm.  Due  to  machining  variations,  however, 
the  as-machined  profiles  differ  from  the  prescribed  profile  in  the  part  drawing  and  were 
measured  using  a  contacting  profilometer  prior  to  being  tested.  The  contact  profilometer  had  a 
vertical  resolution  of  approximately  10  nm.  Two  examples  of  measured  profiles  are  plotted  in 
Figure  2.  Note  that  the  scale  on  the  y-axis  is  significantly  magnified  to  highlight  variations  in  the 
measured  profiles.  It  is  clear  from  these  measurements  that  the  central  flat  section  of  the  pad  is 
not  truly  flat  and  that  the  radii  at  the  edges  do  not  always  match  the  prescribed  radii  either  in 
position  or  sharpness.  Additionally,  there  is  measurable  variability  in  the  profile  geometry  from 
pad  to  pad.  It  was  desired  to  capture  this  variability  through  a  statistical  analysis  so  it  could  be 
used  in  the  subsequent  probabilistic  analysis.  To  achieve  that  objective,  seventy-seven  contact 
pads  were  measured. 

In  order  to  model  the  contact  profile  variability  in  the  probabilistic  analysis,  a 
mathematical  description  of  the  contact  profile  was  needed  that  could  capture  the  important 
features  of  the  profile.  Key  features  were  believed  to  be  the  actual  radii  of  the  edges  of  the  pad 
as  these  can  significantly  affect  the  stress  at  the  edge  of  contact.  Also  important  was  the  flatness 
of  the  central  portion  of  the  pad.  A  rounded  pad  versus  a  very  flat  pad  could  also  significantly 
affect  the  contact  stresses.  A  piecewise  curve  with  twenty-one  parameters  was  defined  to 
represent  the  measured  contact  profiles  and  was  written  as 
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(5) 


rkLx  +  bL  x  <  at 

To  +  V(^L)2  -  0-*o)2  ax<x  <a2 

y  —  \  c0  +  C\X  +  c2x2  +  c3x3  +  c 4x4  +  c5x5  +  c6x6  a2  <  x  <  a3 

yR  +  J(Rr)2  ~  (x-x£)2  a3  <  x  <  a4 

JcRx  +  bR  x  >  a4 

where  x  and  y  represent  a  Cartesian  coordinate  system  with  x  along  the  profile  and  y  is 
the  height  of  the  profile.  In  the  current  problem  the  pad  profile  is  pressed  into  a  flate  plate  so  the 
function  y  is  also  the  gap  function  h.  The  central  segment  of  the  profile  {a2<x<a3)  was 
represented  by  a  sixth  order  polynomial  to  allow  the  nominally  flat  section  to  have  curvature. 

Just  outside  the  central  section  were  two  circular  arcs  to  represent  the  radii  at  the  edges  of 
contact.  Just  beyond  the  radii  were  flat  taper  sections.  Through  various  constraints  at  x  =  a2  - 
a4,  the  number  of  parameters  was  reduced  from  twenty-one  to  thirteen.  These  thirteen  variables 
kL,  kR,  bL.  bR,  yd\  xd\  xR ,  Rl,  Rr,  c/,  c2,  C3,  C4,  c5,  and  c6  were  determined  from  nonlinear 
regression  of  each  of  the  measured  seventy-seven  profiles.  These  variables  were  referred  to  as  ft, 
where  i  ranged  from  0  to  12.  The  remaining  8  variables  in  Equation  1  were  functions  of  y,  as 
defined  by  the  constraints.  The  correlation  coefficients  for  y  are  listed  in  Table  2,  showing  high 
correlation  among  some  of  the  variables.  Finally,  the  results  of  the  nonlinear  regression  from  all 
77  profiles  were  fit  to  a  multivariate  normal  distribution.  The  mean  and  standard  deviation  of  the 
regression  parameters  are  listed  in  Table  1.  These  variable  values  will  result  in  a  pad  profile  or 
gap  function,  h,  with  units  of  pm.  Three  sample  realizations  of  these  fitted  profiles  are  plotted  in 
Figure  6  along  with  the  prescribed  profile.  The  measured  and  the  sampled  profiles  often  differ 
significantly  from  the  prescribed  profile,  particularly  near  the  edge  of  contact. 
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(rarl)  (mrl)  6 


x  (jum) 


Figure  5:  Prescribed  and  sample  measured  pad  profiles. 
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Table  1:  Summary  of  random  variable  statistics. 


Random  Variable 

No. 

Mean 

St.  Dev. 

Distribution  Type 

Initial  crack  size,  a,  (pm) 

20 

15.1 

8.48 

Lognormal 

Friction  coefficient,  juave 

x2 

0.302 

0.021 

Correlated  normal 

Partial  slip  slope,  a 

x3 

1.96 

0.120 

p23  =  -  0.375 

Crack  growth,  logio(Cy) 

X4 

-14.6 

0.486 

Correlated  normal 

Crack  growth,  n  / 

x5 

7.19 

0.715 

p45  —  -0.9973 

Crack  growth,  logi0(C2) 

X6 

-11.8 

0.157 

Correlated  normal 

Crack  growth,  n2 

X7 

3.81 

0.146 

p67  =  -0.9751 

Profile,  y0 

=  kL 

X8 

0.181 

5. 84x1 0'3 

Correlated  normal 

Profile,  /j 

L 

=  yo 

X9 

-2335 

410 

(see  Table  XX) 

Profile,  y2 

=  rl 

Xw 

2333 

411 

Profile,^ 

L 

=  x0 

Xu 

-1612 

37.7 

Profile,  % 

=  Rr 

x12 

2289 

379 

Profile,  y5 

=  X0R 

X13 

1620 

35.4 

Profile,^ 

=kR 

X14 

-0.183 

4.96xl0'3 

Profile,  y7 

=  Cl 

X15 

-2.00xl0'4 

6. 20x1  O’4 

Profile,^ 

=  C2 

X16 

-1.21xl0'6 

l.OlxlO’6 

Profile,  yg 

=  C3 

X]7 

1.53xl0"10 

6.38xl0"10 

Profile,  yio 

=  C4 

Xi8 

9.80xl0"13 

6.16xl0'13 

Profile,  Yu 

=  C5 

Xi9 

-3.77xl0"17 

1.87xl0"16 

Profile,  y12 

=  C6 

X2o 

-3.80xl0~19 

1.59xl0"19 

Table  2:  Pad  profile  regression  correlation  coefficients 


7o 

7i 

72 

73 

74 

75 

76 

7? 

78 

79 

7io 

7n 

7n 

To 

1 

-0.079 

0.078 

-0.119 

-0.051 

0.184 

-0.364 

0.153 

0.175 

-0.107 

-0.130 

0.073 

0.112 

n 

-0.079 

1 

-1.000 

-0.921 

-0.404 

0.321 

-0.122 

0.093 

-0.205 

-0.203 

0.281 

0.275 

-0.246 

72 

0.078 

-1.000 

1 

0.921 

0.404 

-0.321 

0.121 

-0.092 

0.207 

0.203 

-0.282 

-0.275 

0.247 

73 

-0.119 

-0.921 

0.921 

1 

0.325 

-0.372 

0.270 

-0.019 

0.104 

0.136 

-0.201 

-0.247 

0.219 

74 

-0.051 

-0.404 

0.404 

0.325 

1 

-0.899 

-0.204 

-0.102 

0.079 

0.065 

-0.116 

0.002 

0.123 

75 

0.184 

0.321 

-0.321 

-0.372 

-0.899 

1 

-0.033 

0.010 

-0.004 

0.062 

0.021 

-0.147 

-0.061 

76 

-0.364 

-0.122 

0.121 

0.270 

-0.204 

-0.033 

1 

-0.095 

-0.242 

-0.017 

0.142 

0.078 

-0.133 

77 

0.153 

0.093 

-0.092 

-0.019 

-0.102 

0.010 

-0.095 

1 

0.140 

-0.876 

-0.027 

0.627 

0.059 

78 

0.175 

-0.205 

0.207 

0.104 

0.079 

-0.004 

-0.242 

0.140 

1 

-0.099 

-0.869 

0.057 

0.765 

79 

-0.107 

-0.203 

0.203 

0.136 

0.065 

0.062 

-0.017 

-0.876 

-0.099 

1 

0.011 

-0.915 

-0.027 

710 

-0.130 

0.281 

-0.282 

-0.201 

-0.116 

0.021 

0.142 

-0.027 

-0.869 

0.011 

1 

0.027 

-0.944 

7n 

0.073 

0.275 

-0.275 

-0.247 

0.002 

-0.147 

0.078 

0.627 

0.057 

-0.915 

0.027 

1 

-0.031 

712 

0.112 

-0.246 

0.247 

0.219 

0.123 

-0.061 

-0.133 

0.059 

0.765 

-0.027 

-0.944 

-0.031 

1 
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4.  Probabilistic  Analysis 

Monte  Carlo  sampling  was  carried  out  to  determine  the  moments  (mean  and  standard 
deviation)  of  the  cycles-to-failure  distribution.  This  procedure  involved  repeated  generation  of 
realizations  of  the  random  variables  and  execution  of  the  deterministic  fretting  fatigue  algorithm 
to  determine  cycles-to-failure.  The  deterministic  fretting  fatigue  algorithm  was  simplified  as 
much  as  possible  to  minimize  the  computational  time  needed  for  each  set  of  random  variable 
samples.  All  necessary  calculations  were  performed  using  MATLAB  since  both  CAPRI  and  the 
crack  growth  calculation  code  were  developed  using  MATLAB.  This  included  calculation  of  the 
contact  forces  using  the  method  demonstrated  by  Gean  and  Farris  [6],  rather  than  modeling  the 
contact  in  a  FEM  analysis.  The  CAPRI  run  time  was  minimized  by  reducing  the  number  of  Fast 
Fourier  Transform  terms  to  the  minimum  required  and  optimizing  the  solver  parameters  for  this 
problem.  Finally,  a  linear  fit  of  the  bulk  stresses  as  a  function  of  the  contact  forces  was  created 
from  a  series  of  FEM  analyses,  rather  than  using  a  new  FEM  analysis  for  each  Monte  Carlo  run. 
Eliminating  the  FEM  analyses  from  the  calculation  of  fretting  fatigue  life  reduced  the 
computational  time  to  approximately  1.2  s  per  simulation  when  mnning  in  parallel  with  4 
processors  on  an  Intel  Xeon  quad  core  workstation.  This  made  mnning  the  Monte  Carlo  analysis 
with  a  significant  number  of  samples  possible.  The  ensemble  of  cycles  to  failure,  Nf,  results 
were  then  analyzed  to  determine  the  mean  and  standard  deviation.  Sufficient  samples  were 
executed  to  ensure  high  confidence  in  the  computed  moments. 

Probabilistic  sensitivities  play  an  important  role  in  determining  insight  into  the  dominant 
factors  governing  a  probabilistic  analysis  and  provide  information  as  to  potentially  effective 
methods  to  improve  the  reliability.  There  are  a  number  of  methods  in  the  literature  such  as 
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scatter  plots  [20],  Pearson  or  Spearman  correlation  [21],  regression  methods  [22],  and  the  Score 
function  method  [23],  among  others,  that  can  provide  useful  information.  The  following  is  a 
description  of  several  methods  used  in  the  current  work. 

Scatter  plots 

Scatter  plots  are  a  two-dimensional  point  plot  of  the  sample  points  versus  the 
corresponding  response  points.  The  sample  realizations  for  each  random  variable  (X)  and  the 
corresponding  results  (T)  are  plotted  on  separate  axes.  If  a  random  variable  is  not  important,  no 
pattern  should  be  discernable,  that  is,  the  samples  for  X  should  mimic  its  marginal  distribution. 
Conversely,  if  a  random  variable  is  important,  the  pattern  of  realizations  for  X  will  be  distinctly 
non-random.  Scatterplots  are  an  inexpensive  but  qualitative  method. 

Regression  and  correlation 

Linear  regression  (LR)  methods  are  well  known  and  widely  available  tools  for  assessing 
variance  contribution.  LR  approximates  the  relationship  between  response  and  the  random 
variables  as 

y(X)  =  Po  +  'LPinxi)  (6) 

where /is  an  arbitrary  function  of  the  random  variable  Xt.  In  general,  LR  may  contain 
quadratic  and  interaction  terms  of  any  number  of  the  variables,  Xj. 

The  results  from  an  LR  model  provide  an  easy  mechanism  to  estimate  the  sensitivities  of 
the  response  (denoted  Z)  standard  deviation  (oz)  to  the  standard  deviation  of  the  input  parameters 
e.g.,  d(j/jd(jj  ,  or  of  the  response  mean  (pz)  to  the  mean  of  the  input  parameters  (pi),  e.g., 

d/uz  / d jui .  This  sensitivity  in  its  nondimensionalized  form  doz  / dai  *  cri  /  az  or 

cjuz  /  qu,  *  //,  /  pz  provides  an  indication  as  to  the  relative  importance  of  the  random  variables 

and  can  be  used  to  estimate  design  improvements. 
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One  of  the  simplest  linear  regression  models  in  terms  of  the  original  variables  is 


Z  —  Aq  +  +  •••  +  AnXn 


In  order  to  simplify  the  sensitivity  equations  derived  below,  Equation  1  can  be  rewritten 


as  the  standardized  model 


z~z  (zi  Zlhl  i  |  X  ~M„) 

n  1  n 


where  Z  denotes  the  sampling  mean  of  the  response  and  Sz  represents  the  sampling 


standard  deviation  of  Z.  Equation  2  can  then  be  rewritten 


z  —  A^Xi  +  •••  +  AnXn 


where  the  conversion  between  the  original  and  standardized  coefficients  is 


A  =^-A. 


Aq  —  Z  —  A1fl1 


AjiMn 


It  is  known  that  for  a  linear  model  of  the  form  Z  =  A0  +  A1X1  +  — h  AnXn  that  the 
expected  value  or  mean  value  of  Z  can  be  written 


E(z)=YJAiE(X,) 


and  the  sensitivities  or  partial  derivatives  are  simply 


The  variance  of  Z  can  be  written  as 


Var(z)=YYjAiAjPij(jicjj 


[24]  and  the  sensitivities  or  partial  derivatives  are  [25] 
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(15) 


^=^-±A 


do:  <7 


jPijVj 


Z  7  =  1 


For  independent  variables  (py  =  0)  this  reduces  to 


dcJz  _  Oj  A2 
do,  o7 


(16) 


Note,  based  on  Eq.  (10),  the  sensitivity  must  be  positive,  whereas,  based  upon  Eq.  (9),  the 
sensitivity  can  be  negative  for  negative  correlation  coefficients. 

Rewriting  the  sensitivity  in  terms  of  standardized  regression  coefficients  yields 


ooi  cr  oz  j=l  Oj  cr,  j=l 


(17) 


Nondimensionalizing  yields 

°  i  _  °z 


S,.= 


do,  o7  o7  o .  1  " 

l  Z  Z  /  /=! 


=  A^AjPij 


(18) 


7=1 


For  independent  variables  this  reduces  to 


s  = 

do ,  O  y 


(19) 


Finite  Difference 

The  finite  difference  (FD)  method  [26]  was  used  in  this  study  to  compute  the 
probabilistic  sensitivities  as  a  method  of  verifying  other  much  more  efficient  methods  such  as 
LR.  FD  is  conducted  by  simply  applying  a  small  change  in  a  random  variable  parameter  (mean 
or  standard  deviation)  and  then  re-running  the  Monte-Carlo  analysis  to  obtain  a  new  response 
distribution.  The  change  in  the  modes  of  the  response  distribution  divided  by  the  change  in  the 
random  variable  parameter,  i.e.  Acrz/Acri,  gives  an  estimate  of  the  probabilistic  sensitivity, 
doz/do^ .  Calculation  of  the  probabilistic  sensitivities  using  FD  can  be  very  computationally 
costly  due  to  the  large  number  of  Monte  Carlo  simulations  required.  In  the  current  problem  there 
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are  20  random  variables  each  with  2  parameters  (mean  and  standard  deviation)  multiplied  by  a 
sample  size  of  50,000  for  each  parameter  results  in  2  million  simulations  or  approximately  1 
month  computation  time.  50,000  was  chosen  as  a  necessary  balance  between  accuracy  and 
computational  cost.  Unfortunately  it  is  unlikely  to  be  enough  samples  to  achieve  convergence  in 
the  sensitivity  estimate.  In  fact,  additional  simulations  were  performed  to  double  the  sample  size 
to  100,000  for  the  FD  analysis  of  several  variables  of  interest  proving  the  estimates  were  either 
converged  (less  than  5%  change)  or  in  some  cases  not  converged.  However,  of  those  variables 
of  interest  that  were  not  converged  changes  in  the  estimates  with  sample  size  of  100,000  were 
modest,  i.e.  30%  or  less.  A  percent  change  in  sensitivity  estimates  for  many  variable  parameters 
were  not  considered  since  the  estimates  were  nearly  zero. 

5.  Results  and  Discussion 
Probabilistic  Life  Prediction 

The  probabilistic  fretting  fatigue  life  prediction  tool  was  exercised  at  several 
experimental  load  cases.  Results  were  generated  at  the  5  applied  loading  conditions  (Fmax  =18, 
20,  22,  24,  and  30  kN)  used  in  the  dovetail  testing  program.  Lives  were  generated  at  each 
loading  condition  initially  using  a  sample  size  of  10,000  for  the  random  variables  initial  crack 
size,  coefficient  of  friction,  partial  slip  slope,  crack  growth  curve,  and  the  pad  profile.  These 
results  are  plotted  on  lognormal  probability  paper  in  Figure  7  along  with  the  experimental  lives. 
The  solid  lines  are  the  distribution  of  the  Monte  Carlo  predicted  failure  lives  with  the  applied 
highest  load  (30  kN)  on  the  left  and  the  lowest  applied  load  (18  kN)  on  the  right.  The 
experimental  lives  shown  were  as  follows:  1,371,000,  10,000,000,  and  693,000  cycles  at  18  kN; 
995,000  and  919,900  cycles  at  20  kN;  249,900,  346,200,  and  497,800  cycles  at  22  kN;  164,000 
cycles  at  24  kN;  and  105,000  cycles  at  30  kN.  These  data  were  quite  limited  for  an  evaluation  of 
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a  probabilistic  life  prediction,  however,  they  were  still  useful  to  ensure  the  predictions  were  near 
the  test  results.  The  comparison  in  Figure  7  shows  that  at  the  higher  loads  the  predictions  were 
closer  than  at  the  lower  loads.  Not  surprisingly,  all  of  the  experimental  results  fall  within  the 
scatter  of  the  simulated  results.  At  the  lower  applied  loads,  however,  the  test  results  seem  to 
deviate  from  the  predictions  by  having  longer  lives.  In  fact,  the  longest  life  at  18  kN  is  over  10 
million  cycles  (test  was  a  run-out)  which  was  at  the  99%  probability  level  of  the  prediction.  If 
the  three  18  kN  experiments  do  belong  to  the  simulated  failure  life  distribution  there  is  less  than 
a  3%  chance  that  this  long  life  would  occur  in  testing.  This  shows  that  the  fracture  mechanics 
life  prediction  model  is  likely  missing  some  of  the  physics  of  the  problem.  It  supports  the  view 
that  crack  nucleation  or  formation  may  be  a  significant  portion  of  the  life  at  lower  stress  levels, 
even  for  a  crack  propagation  model  has  an  initial  crack  size  at  a  microstructural  scale. 
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Crack  Growth  Life,  N  (cycles) 

Figure  7:  Probability  plot  of  experimental  (points)  and  predicted  (lines)  lives  at  several  applied 

loads. 

The  results  plotted  in  Figure  7  were  generated  with  just  10,000  samples  to  minimize 
processing  time  for  each  condition  analyzed.  Since  the  number  of  samples  was  relatively  small, 
an  understanding  of  the  variance  in  the  estimates  of  mean  and  standard  deviation  of  Nf  was 
important.  Variance  estimates  of  the  mean  and  standard  deviation  of  the  distribution  were 
calculated  at  each  of  the  applied  forces  and  listed  in  Table  3.  Since  the  distribution  of  Nf  was 
fairly  lognormal,  the  modes  of  the  distribution  of  the  logarithm  of  Nf  were  calculated.  The 
estimates  for  mean  of  the  logio  of  Nf  were  then  converted  back  to  cycles.  The  columns  showing 
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the  95%  confidence  bounds  of  the  mean  were  then  converted  back  to  cycles  or  the  ratio  of  the 
lower  bound  (LB)  and  upper  bound  (UB)  to  the  mean  cycle  count.  The  coefficient  of  variance 
(COV)  was  also  calculated  entirely  in  log  units.  Interestingly,  the  COV  decreased  with  increased 
applied  loading,  which  was  consistent  with  fatigue  test  results  in  Ti-6A1-4V  [18]  and  it  is 
generally  true  in  metallic  aerospace  materials  that  at  higher  stresses  the  variance  in  fatigue  life 
decreases.  Also,  in  Table  3  the  estimates  of  the  mean  and  standard  deviation  of  the  distribution 
were  calculated  for  smaller  and  larger  numbers  of  samples  at  Fmax  =  22  kN.  With  more  samples, 
the  confidence  bounds  shrink  significantly  as  expected.  In  addition  to  quantifying  the  modes  of 
the  distribution,  PDF’s  were  fit  to  the  distribution.  Figure  8  is  a  plot  of  two  PDF  types  fit  to  the 
Nf  results  at  Fmax  =  22  kN.  The  solid  line  is  a  nonparametric  kernel  estimate  of  the  PDF  [26], 
while  the  dashed  curve  was  a  lognormal  fit  to  the  Nf  distribution.  The  differences  in  the  PDF’s 
showed  that  the  distribution  of  failure  lives  was  skewed  toward  longer  lives  and  also  had  a  higher 
peak.  This  skewness  toward  longer  lives  in  the  tail  can  also  be  observed  in  Figure  7,  and  this 
trend  increased  with  smaller  applied  loads. 
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105  106 
N  (cycles) 

Figure  8:  PDF  estimates  of  N/ for  10,000  simulations  at  Fmax  -  22  kN. 


Table  3:  Mean  and  coefficient  of  variation  (c.o.v.)  of  the  Monte  Carlo  analysis  results  including 
the  95%  upper  bound  (UB)  and  lower  bound  (LB)  confidence  limits 


Mean  of  logio (Nf)  COV  of  log10(/V/) 


Fmax  (kN) 

Sample 

Size 

LB  /  Mean 

Mean 

(cycles) 

UB  /  Mean 

LB / COV 

COV 

UB / COV 

18  kN 

10,000 

0.989 

872,600 

1.012 

0.986 

4.26% 

1.012 

20  kN 

10,000 

0.990 

509,900 

1.010 

0.988 

4.01% 

1.012 

22  kN 

1000 

0.973 

328,400 

1.028 

0.960 

3.54% 

1.042 

22  kN 

10,000 

0.991 

326,200 

1.009 

0.987 

3.72% 

1.013 

22  kN 

500,000 

0.999 

324,300 

1.001 

0.998 

3.70% 

1.002 

24  kN 

10,000 

0.992 

225,400 

1.008 

0.989 

3.51% 

1.014 

30  kN 

10,000 

0.993 

95,400 

1.007 

0.988 

3.24% 

1.012 

Sensitivity  Results 

Perhaps  the  simplest  method  to  qualitatively  determine  sensitivity  is  through  a  scatter 
plot.  Figure  9  below  shows  the  scatter  plots  that  relate  cycles-to-failure  ( y  axis)  with  each 
random  variable  (x  axis)  for  10,000  samples.  The  scatter  plot  for  random  variable  X2  shows  a 
definite  pattern.  Correlation  coefficients  (Pearson  or  Spearman)  for  each  variable  relate  to  the 
amount  of  variance  in  Y  that  can  be  apportioned  to  X,.  The  results  for  Pearson  correlation 
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coefficients,  given  in  Table  4,  indicate  that  variable  X2  (coefficient  of  friction)  is  by  far  the 
dominant  variable,  contributing  65%  of  the  variance  in  Nf.  It  must  be  pointed  out,  however,  that 
the  correlation  coefficients  are  obtained  without  consideration  for  the  simultaneous  variations  in 
other  random  variables.  A  better  estimate  of  variable  importance  can  be  determined  using  linear 
regression,  discussed  below. 
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Figure  9:  Scatter  plots  of  the  model  output  cycles  to  failure,  Nf,  versus  the  input  random 

variables,  Xj  through  X2o 

In  this  research,  an  LR  model  of  the  form  of  Equation  6  was  used,  where  X  denotes  a 
vector  of  random  variables,  y  represents  the  cycles-to-failure  and  the  /?’ s  are  coefficients  that  are 
fit  to  the  analytical  results.  Table  5  shows  the  results  of  a  “best  model”  fit  using  a  specified 
number  of  variables.  For  example,  using  only  2  variables,  the  combination  of  X2  and  Xk,  account 
for  the  largest  percentage  of  the  variance  in  Nf.  Note,  once  X2  is  included  in  the  model,  adding 
Xi6  adds  17  percent  to  the  R"  sum,  however,  the  R"  value  using  the  Pearson  correlation 
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coefficient  for  X2  without  considering  any  other  random  variables  is  25%.  From  the  results,  one 


can  see  quickly  the  diminished  returns  offered  after  the  first  few  random  variables.  For  example, 
if  only  5  random  variables  are  used  in  an  LR  model,  this  model  would  account  for  85%  out  of  a 
possible  91%  of  the  output  variance.  It  is  somewhat  surprising  that  a  simple  linear  model  with 
respect  to  the  random  variables  can  account  for  such  a  large  percentage  of  the  output  variance. 

Table  4:  R~  values  for  the  random  variables 


Variable 

R~ 

20 

0.024 

*2 

0.645 

V, 

0.093 

20 

0.001 

Xs 

0.002 

x6 

0.006 

X7 

0.002 

20 

0.002 

Xg 

0.004 

Xjo 

0.004 

Xn 

0.003 

X,2 

0.005 

X13 

0.017 

X,4 

0.008 

X,5 

0.007 

Xl6 

0.063 

X17 

0.013 

Xj8 

0.020 

X19 

0.020 

X20 

0.038 
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_ Table  5:  Best  model  linear  regression  results 

~N  IF  C(p)  Random  Variables  in  Model 


1 

0.645 

28,370 

z2 

2 

0.713 

21,020 

X2 

X16 

3 

0.754 

16,620 

Z2 

X9 

Xjo 

4 

0.824 

9000 

X2 

X16 

X18 

X20 

5 

0.846 

6660 

X; 

X2 

X16 

X18 

X20 

6 

0.855 

5710 

X, 

X2 

Xj3 

X16 

X18 

X20 

7 

0.874 

3640 

Xj 

X2 

X6 

X7 

X]6 

X]8 

X20 

8 

0.883 

2690 

X , 

X2 

X6 

X7 

X,3 

X,6 

X18 

X20 

9 

0.898 

1060 

Xj 

X2 

X6 

X7 

X16 

X]7 

X]8 

Xjg 

10 

0.901 

720 

X , 

X2 

X6 

X7 

x13 

X,6 

Xl7 

X18 

Table  6  shows  the  results  using  LR  considering  natural  groupings  of  random  variables. 
The  variables  were  partitioned  into  4  groups:  (X/  -  initial  crack  size,  X2-X3  coefficient  of  friction 
-  partial  slip  slope,  X4-X7  crack  growth  parameters,  X8-X20  geometry  profile).  The  results 
indicate  that  group  2  is  dominant,  followed  by  the  geometry  profile.  Surprisingly,  the 
traditionally  dominant  random  variables  in  fatigue,  crack  growth  rate  and  initial  crack  size  are 
not  significant  in  terms  of  their  contribution  to  the  output  variance,  e.  g.,  these  parameters  could 
be  modeled  deterministically. 

Table  6:  Linear  regression  group  analysis  showing  the  effect  of  each  random  variable  group  on 

the  model  R2 


Group 

RV  Added 

#  RV 

Model 

R2 

R1 

C(p) 

F  value 

1 

x2,x3 

2 

0.645 

0.645 

28,370 

9080 

2 

X8-X20 

15 

0.210 

0.855 

5710 

1110 

3 

x4-x7 

19 

0.030 

0.885 

2420 

660 

4 

X , 

20 

0.022 

0.908 

20 

2410 

Next,  we  wish  to  quantify  the  probabilistic  sensitivities,  particularly  the  sensitivity  of  the 
standard  deviation  of  the  response  (life)  to  the  random  variable  parameters  mean,  daz/dph  and 
standard  deviation,  daz/da t.  Figure  10  is  a  bar  chart  of  the  normalized  sensitivity  of  the 
standard  deviation  of  the  cycles  to  failure,  Nf,  response,  denoted  Z,  to  the  standard  deviation  of 
the  random  variables,  Xj,  simply  denoted  by  i  in  Equation  19.  Results  from  the  finite  difference 
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analysis  and  linear  regression  show  very  good  agreement.  The  analysis  shows  that  the  standard 
deviation  of  the  response  is  by  far  most  sensitive  to  the  standard  deviation  of  the  coefficient  of 
friction  (X2)  and  to  several  parameters  of  the  flat  section  of  the  contact  profile  (AV,  -  X20).  The 
slope  of  the  first  segment  of  the  bilinear  crack  growth  law  (X5)  and  the  intercept  of  the  second 
segment  (Ag)  also  have  relatively  high  sensitivities.  Also  notable  is  the  sensitivity  of  the  initial 
crack  size  (A/)  is  quite  small.  The  sensitivities  of  these  variables  is  essential  to  understand  when 
evaluating  which  variables  are  most  important  to  the  process  being  analyzed.  One  must  also 
consider,  however,  the  likelihood  that  the  estimates  of  these  random  variable  parameters  will 
change  due  to  additional  or  better  quality  data,  or  that  these  parameters  can  change  through  an 
alteration  to  the  system,  i.e.  machining  tolerances  or  addition  of  coatings. 


Random  Variable,  X, 


Figure  10:  Comparison  of  5,  for  FD  and  LR. 


Figure  10  helped  to  identify  that  the  standard  deviation  of  A;<j,  X/g,  and  X20  were 
parameters  with  a  relatively  high  sensitivity.  It  is  not  intuitive,  however,  how  these  variables 
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affect  the  pad  profile.  Further  analysis  of  the  profile  has  shown  that  the  flatness  of  the  “flat” 
section  ( a2  <  x  <  a  3)  of  the  pad  profile  is  strongly  affected  by  the  values  of  these  variables.  Here, 
flatness  is  defined  as  the  peak  to  peak  height  difference  in  y.  A  probabilistic  fretting  analysis 
was  performed  the  pad  profile  being  the  only  group  of  random  variables.  The  remaining 
variables  were  deterministic.  The  results  showed  a  significant  correlation  between  log10  Nf  and 
X/6,  X/8,  and  X20  of  0.35  to  0.55,  however,  the  correlation  between  log10  Nf  and  flatness  was 
0.77.  Figure  11  is  a  scatter  plot  showing  this  relationship. 


max  min 

Figure  11:  Scatter  plot  of  fatigue  crack  growth  life  versus  profile  flatness  showing  far  stronger 
correlation  than  with  any  of  the  random  variables. 


6.  Conclusions 

A  new  probabilistic  fretting  analysis  was  developed  to  investigate  the  relative  importance 


of  typical  fretting  input  variables  on  the  predicted  failure  lives.  Several  random  variable  inputs 
were  identified  and  PDF’s  were  quantified  using  laboratory  data.  Monte  Carlo  sampling  of  the 
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input  PDF’s  was  performed  and  a  deterministic  analysis  was  repeatedly  run  using  the  sampled 
inputs  to  obtain  a  distribution  of  predicted  fretting  lives.  The  results  were  compared  to 
experimental  fretting  fatigue  tests  and  the  predictions  were  much  closer  at  higher  stresses  than  at 
lower  stress  possibly  due  to  substantial  crack  nucleation  lives  at  the  lower  stresses.  The  results 
showed  that  considerable  scatter  in  fretting  lives  can  be  expected  based  on  variability  in  the 
material  properties,  contact  profiles,  coefficient  of  friction,  and  contact  force  response. 
Interestingly,  the  dominant  variables  in  terms  of  contribution  to  the  cycles-to-failure  variance 
were  the  coefficient  of  friction  and  several  terms  within  the  geometry  profile;  whereas,  the 
traditionally  dominant  variables  in  fatigue  crack  growth  analyses,  initial  crack  size  and  crack 
growth,  were  not  as  significant. 
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